from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm


plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题
plt.rc('font', size=10)
plt.rc('font', family='SimHei')

m_A = 4866
m_B = 2433
f = 4890
delta_m = 1165.992
omiga = 2.2143
k_3 = 167.8395  # 兴波
k_5 = 80000  # 弹簧劲度系数
k_6 = 10000  # 阻尼系数
ro = 1025
g = 9.8
s = np.pi

def P(k_6):
    # 定义一个方程组（微分方程组）
    def pfun(y, t):
        y0, y1, y2, y3 = y
        return np.array([
            1 / (m_A + delta_m) * (
                        -k_3 * y0 - ro * g * s * y1 + f * np.cos(omiga * t) + k_5 * (y3 - y1) + k_6 * (y2 - y0)),
            y0,
            1 / m_B * (-k_5 * (y3 - y1) - k_6 * (y2 - y0)),
            y2
        ])

    t = np.arange(0, 500, 0.2)  # 创建自变量序列
    soli = odeint(pfun, [0, 0, 0, 0], t)

    re = 0
    for i in range(200, len(soli[:, 2])):
        re += k_6 * (soli[:, 0][i] - soli[:, 2][i]) ** 2 * 0.2
    return re / 460

y = []
k = [i for i in range(1, 100100, 100)]
for i in tqdm(range(1, 100100, 100)):
    y.append(P(i))
print(y)
plt.figure(figsize=(6, 4))
plt.title('平均输出功率与阻尼系数曲线')
plt.plot(k, y, label='P(k)', linewidth=3.0)
plt.xlabel('阻尼系数')
plt.ylabel('平均输出功率')
plt.legend()
plt.grid(True, which="both", ls="--")
plt.show()
